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We analyze the effect of tlie long-range interaction on the transport properties through ordered 
and disordered one-dimensional metallic nanoparticle arrays. We discuss how the threshold voltage, 
the I-V curves and the voltage drop through the array are modified as compared to the case in 
which interactions are restricted to charges placed on the same island. We show that some of 
these modifications are due to finite interactions between charges in different nanoparticles while 
other ones are due to interactions between charges in the islands and those at the electrodes, what 
produces a polarization potential drop through the array. We study the screening of the disorder 
potential due to charged impurities trapped in the substrate and find that long-range interactions 
introduce correlations between the disorder potentials of neighboring islands. 

PACS numbers: 73.23.-b,73.63.-b,73.23.Hk 



I. INTRODUCTION 



Nanoparticle arraysi^^>^>&>i>a>a>ia^i^^^^i^^^^> 

are a perfect system to analyze correlated elec- 
tronic transport and have received a lot of exper- 
imental and theoretical attention during the last 
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In spite of this, their properties and I-V characteristics 
are not well understood. The presence of disorder in 
these arrays complicates the analysis^^^i^ . The two 
most relevant and experimentally unavoidable types of 
disorder are charging disorder due to random charges 
quenched at the substrate and resistance disorder due to 
small changes in the distance between the nanoparticles 
and the exponential dependence of the resistance on the 
distance. Voids in the lattice or a non-homogeneous 
distribution of nanoparticle size can be also presenl;^ but 
this source of disorder is less important in some recent 
experiments. 

Theoretical analysis have mainly considered arrays in 
which each nanoparticle is capacitively coupled only to 
its nearest neighbor o^'^'^'^i'^^'^^i^^'^'^i^^ , especially the case 
in which this coupling is small. The truncation of ca- 
pacitive coupling to nearest neighbor results in an inter- 
action between charges in different conductors which de- 
cays exponentially with the distance between then*^!^. 
This limit is relevant for those arrays coupled to a gate 
electrode^, as the mobile charges in this lead effectively 
screen Coulomb interactions. A complete description of 
the non-equilibrium transport through arrays has been 
presented only very recently^ and is restricted to one- 
dimensional metallic nanoparticle arrays in which inter- 
actions are finite only when charges are placed on the 
same conductor. 

Self-assembled arrays fabricated nowadays are de- 
posited onto insulating substrates and generally lack a 
gate voltage. In these arrays, the screening of long-range 
interactions is less effective, but the proximity of other 



conductors, both islands and leads, modifies its value 
compared to a 1/r Coulomb law^i^. The electrodes 
20.EDBi ribute to the screening of the interaction. Theoreti- 



cal analysis including the effect of long-range interactions 
are scarce and limited^^i^i^ to numerical results or par- 
ticular cases. In this paper we perform a detailed analysis 
45 .a£.th e effect of the long-range character of the interaction 



on the transport properties through ordered and disor- 
dered one-dimensional arrays. We consider the influence 
of charging and resistance disorder. The interactions are 
described by an inverse capacitance matrix, in which cou- 
pling to other conductors is not truncated. The matrix is 
calculated including the effect of screening and the cur- 
rent is computed numerically. Some analytical approx- 
imations are also discussed and compared with the nu- 
merical results. This comparison allows us to understand 
the origin of the behavior found. To analyze the trans- 
port we put forward a model to describe the electrostatic 
interactions among the charges occupying the islands and 
the electrodes, which allows for a clear description of the 
relevant contributions to the transport. In particular, 
we introduce the concept of polarization potential drop 
at the junctions between the nanoparticles. As in our 
previous study of the short-range caso^, we analyze the 
threshold voltage, flow of current and voltage drop for 
the cases with and without charge or junction resistance 
disorder. 

The organization of the paper is as follows. Section 
II describes the system under study and the model used 
to analyze the transport properties and compares it with 
previous models used in the literature. In section III 
we discuss the correlations introduced by the long-range 
character of the interaction in the disorder potential. Sec- 
tion IV, V and VI respectively analyze the threshold volt- 
age, I-V characteristics and voltage drop through the ar- 
ray. In section VII we summarize our results. Finally, 
in Appendix I we describe the two methods employed to 
calculate the screened interaction between charges, and 
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in Appendix II we discuss the effect of screening on the 
distance-dependence of this interaction. 



II. THE MODEL 

We consider an array composed of N metaUic spheres 
of radius r'*' and center to center distance 2r**' + d. 
Throughout lengths are measured in units of r"' and 
energies in units of the charging energy of an isolated 
nanoparticle Elf'- = 1/(2C"'''), with C"*' the nanoparticle 
capacitance when it is isolated. Here and in the follow- 
ing the electronic charge e = 1. As d/r'*' decreases, the 
charges in an island feel more the effects of charges in 
other islands. In this sense we say that decreasing the 
spacing d/r**' increases the range of interactions. Ex- 
perimentally arrays in which the metallic nanoparticles 
are capped with thiols have d/r^^' ~ 0.5. Arrays self- 
assembled via other types of molecules, like DNA, allow 
for the investigation of larger d/r*'*' values. 

We model each nanoparticle with a continuum level 
spectrum. The tunneling barriers which separate them 
have a resistance much larger than the quantum of resis- 
tance and we treat the transport at the sequential tun- 
neling level^^. To analyze the transport the array is sand- 
wiched between two large electrodes. We assume that the 
electrodes are not ideal voltage sources, but have finite 
self-capacitances. The potential on the leads will thus 
fluctuate in response to all tunneling events, even those 
that do not directly involve the electrodes. For a more 
extensive discussion of these assumptions see^ 

The energy of the system is given by 
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Labels and N + l refer to source and drain electrodes 
and 1, . . . , iV to the islands. Latin capital and lower case 
letters are used to denote electrodes and islands respec- 
tively. Greek indexes will be used when the labels refer 
to both islands and electrodes, (pf"^" is a random potential 
at each island, present in charge disordered arrays, and 
zero in the clean case. Charges Qq and Qn+i maintain 
source and drain electrodes at potentials Vq and Vat+i, 
respectively. 

The electrostatic interactions in our system are de- 
fined through an inverse capacitance matrix that 
directly includes all the conductors in our system. All 
the elements are positive. The inverse capacitance 
matrix is symmetric, C~p = C*^^, and has dimension 
(A^ + 2) X (A^ + 2), i.e. it includes the electrodes. We have 
developed two numerical methods to calculate the inter- 
action potential (inverse capacitance matrix) of an array 
of spheres. These methods are explained in Appendix 
I. The properties of this interaction, how it compares 
to previous calculations for the case of cubic, cylindrical 
(with array and island axis coUinear) and disk shaped 
islands (with array and island axis perpendicular) and 



the screening effect of the electrodes are discussed in Ap- 
pendix II. For simplicity we have considered two large 
spherical leads. The coupling between the islands and 
the electrodes is included through C~j^ and C^\. 

The current is computed numerically by means of a 
Monte Carlo simulation, described i n^3,56 jg controlled 
by the probability of a tunneling process, given by 



r(A£;) 
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R exp{AE/KBT) - 1 



(2) 



Here R is the resistance of the junction through which 
tunneling takes place and AE is the difference between 
the energy of the system before and after the tunneling 
event. In the following we restrict the discussion to zero 
temperature. Then r(A£:) = -AE/RtQ{-AE). It is 
finite only when AE is negative. 

After a bit of algebra, the change in energy due to a 
tunneling event from a to (3 can be rewritten as 
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Here, the excitonic energy E^ ^ necessary to create an 
electron-hole pair in an uncharged array is given by 
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The second term in (|3]) can be seen as the potential differ- 
ence between the sites involved in the tunneling. At the 
electrodes (j)o = Vo = aV and 4'n+i — Vn+i — {a — 1)V, 
where V is the bias potential and we have introduced the 
bias asymmetry parameter a as in refi^. Both a — 1/2, 
also denoted as symmetric bias, and a = 1 have been 
used in the literature, a characterizes how the bias volt- 
age is partitioned between source and drain chemical po- 
tential shifts. Since no physical properties depend on the 
overall zero of energy, varying a in our model is entirely 
equivalent to rigidly shifting all impurity potentials by 
—aV. Since in our model all transport occurs by trans- 
fer between adjacent nanoparticles, the evolution of a 
nanoparticle array as the bias voltage is applied is sensi- 
tive to a, and we believe that the dependence on a could 
in principle be observable, see discussion in^. 

At the islands the potential can be decomposed into 
three terms = (j>f''^ + 0f °' -I- the disorder potential 
(fif^^ due to random charges in the substrate , the polar- 
ization potential 0^°' at the island induced by the elec- 
trodes at finite bias and the potential due to the charges 
in the nanoparticles 0^^. Here 
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junctions and independence on a. From ([6]) we see that a 
priori A" depends both on the geometry of the electrodes 
and the array and on how this is biased. It depends on 
how the charges in the islands and the electrodes inter- 
act. For most capacitance matrices, in particular for the 
capacitance matrices discussed here, the polarization po- 
tential is not linear in the island label i and the potential 
drops are larger close to the biasing leads. In the on- 
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C^^ can be interpreted as a modification of the interac- 
tion between the charges in the islands due to the prox- 
imity of the electrodes at a fixed potential. For the case 
i — j m. which both charges are on the same island this 
modification was already discussed in^, as the interac- 
tion of a soliton with a passive edge. Expression ([9]) 
shows that not only when the charges are in the same 
island, but also when they occupy different islands, their 
effective interactions are modified by the presence of the 
voltage-biased leads. Two types of terms can be differ- 
entiated in the modification of this interaction. The last 
two terms in Q or direct terms, can be viewed as the 
interaction between a charge in island i and the image 
charge at one of the electrodes induced by the charge 
in island j. This term is affected by the presence of the 
other electrode. On the other hand, the terms containing 
C'o/v-i-i' or indirect terms, reflect the interaction between 
the image charges in both electrodes. Direct and indi- 
rect terms have opposite sign. The direct term reduces 
the effective interaction; the indirect one increases it. We 
emphasize that ([31) to ([9]) follow from ^ after straight- 
forward and trivial algebra. We have just defined a few 
quantities and split the change in energy AE' and poten- 
tial 4>a in several terms to facilitate the physical inter- 
pretation of the transport properties. 

We can also define the potential drop at each junction 



(10) 



with the corresponding disorder, polarization and charg- 
ing terms (j)''"*^ and <I>^''. Label i for a junction runs 
from 1 to + 1 and refers to the one between conduc- 
tors i — \ and i. The polarization potential drop at each 
junction. 
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does not depend on the resistance of the junctions, but 
on the electrostatic interactions of each island with the 
voltage-biased leads. Here Xq = a and A^_|_;^ = a — 1. 
A linear drop of the polarization potential implies 

(/)f°' ^ [a- V, requires = l/(^ + 1) for all 



site case, discussed in ref;5S A 
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is finite only at 



junctions 1 and + 1 and given by A^f°"***'^J — q, ^nd 

A^^^"***'^-'^ = a — 1. In general, as the range of the interac- 
tions between charges increases. A" is more homogeneous 
In previous models^ had dimension N x N and 
the inverse self-capacitances of the electrodes, Cqq and 
C 
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N+1,N + 



I , and the inverse mutual capacitances between 
them, CqIj_^^ = C^^]^ q, were neglected. In our model 
small quantities whose values will not significa- 
tively affect the quantitative results. These other models 
do not include the indirect term in ([9]) either. On the 
other hand, an expression analogous to (O can be de- 
fined in other models. For example, in the model previ- 
ously discussed by Middleton and WingreenSl, with an 
N X N capacitance matrix, the interaction between an 
island and an electrode is given by the interaction of the 
charge of the island with charges induced by the electrode 
in the islands immediately adjacent to the electrode. This 
interaction results in a polarization potential character- 
ized by 
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Here Ci-ei is the capacitance between the source or drain 
electrode and an adjacent island. In this model, except 
in the extreme long-range case, see below, the polar- 
ization potential does not decay linearly with distance 
and depends on the asymmetry of the bias potential. 
Other models, however, impose a uniform polarization 
drop through the array^ 

As discussed in^ and section VI, even if the polariza- 
tion potential does not drop linearly and is independent 
on the junction resistance, the average total potential 
drop depends on the resistance (via the average charge 
occupation of the islands) and for homogeneous resis- 
tances a linear drop is partially recovered at large bias 
voltages. 

Whenever not specified we assume that all the junction 
resistances Ri are equal and given by Rt- The effect 
of non homogeneous resistances will be studied in two 
ways. In the first case, one of the junction resistances 
at a given position is larger than the other ones (given 
by Rt)- In the second case the value of the resistances, 
varying in between two values is randomly assigned to 
the junctions. To mimic that disorder in resistances orig- 
inates in variations in distances between the islands and 
the exponential dependence of the junction resistance in 
the distance between islands the junction resistance is 
given by i? = Roexp{"fdist) with Rq and 7 input pa- 
rameters and dist = 1 + random/ 2. Here random is a 
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random number between and 1. In the paper, we have 
used Ro = 1.1825i?r and 7 = 1.526,1.95,2.84. With 
these values the resistance changes respectively between 
(5-ll)i?T, (8-21)i?T and (23-83)i?T. 



III. SCREENING OF DISORDER POTENTIAL 

Charged impurities trapped in the substrate under- 
lying the nanoparticle array create random potentials 
at the nanoparticles. In molecularly assembled arrays, 
charge transfer to the organic molecules surrounding the 
nanoparticles results in non-integer random charges at 
the islandg^*^ . Charging disorder is included in our model 
through a random potential at each island In prin- 
ciple, (j)f^'' take values larger than the charging energy 
El^K However, for large values of the disorder potential, 
charges flow to compensate for these large fluctuations. 
In this section we analyze the effect of the long-range in- 
teraction on the final distribution of disorder potential 
as compared to the case with onsite interactions. We 
find that the distribution of probabilities P{(j>'^''^) and 
pj-(j)dis-j modified. The maximum and minimum val- 
ues of and {^f^^} are modified compared to the 
short range case and given by ±C~^^ /2 and E^^^ . Cor- 
relations between the disorder potentials of neighboring 
islands are introduced. 

If interactions between the charges are short-range, 
(Cj^^ = 5ij), the set of disorder potentials once 
the screening of the potential due to the mobile charges is 
taken into account, is uniformly distributed in the inter- 
val -Ei''' < (jif"" < El'K The probability associated with 
each pair, (c^^**, , is a constant, see Fig. 1 and the 

distribution of the probabilities of the potential drops due 
to disorder across the array j mictions, = — (^fl'^ , 
has the form^ 

1-4^ (13) 

and — iEl"'-. In the presence of long-range in- 

teractions, the charges which flow to compensate the 
large fluctuations of the disorder potential, influence the 
value of the total potential at neighboring islands. As a 
consequence, the screened disorder is correlated^. The 
probability of each pair {(j)f^^ , (j)f^i) is no longer a con- 
stant. P($''*'') depends on the inverse capacitance matrix 
C^^. In order to analyze these correlations and obtain 
the proper disorder potential distribution we assign the 
potentials by first randomly assigning potentials to the 
islands (t)f^-b^^^^ in the interval -W < cf)f < W 
with W larger than the charging energy. We then find 
the equilibrium configuration of charges {Qj'^} that oc- 
cupy the array with island disorder potentials 
and grounded leads {Vq = Va^+i = 0) and redefine the 



potentials at each site using the expression 

N 

The effect of the screening charges {Qj"^} is included in 
the redefined potentials {^f*^} so we then reset the num- 
ber of charges at each site to zero to avoid doublecount- 
ing the charge when we calculate the total electrostatic 
energy of our system. 

Following the redefinition of the disorder potentials, 
we find that on average the distribution of the disorder 
potentials {<j)f^^} and the disorder potential drops {<&f"} 
between adjacent islands are independent of W. The 
values of {0f and {$f^} are bound by ±C-,^/2 and 
±Kf respectively, as the total energy of the system is 
at a global minimum when the original disorder configu- 
ration ^(j)d-ts-bare^ Screened out by the charges {Qj"} ■ 
When in this state, adding an additional charge to any is- 
land in the array increases the energy of the system. The 
energy of adding an additional charge to an island from 
a large electrode outside the system with neglible self in- 
verse capacitance is given by Ef'^'^ = {l/2)C~l^ ± (j)f^'^ 
where the top (bottom) sign refers to the change in en- 
ergy of the system as a result of adding a positive (nega- 
tive) charge to island i. Since 

^add > when the array is 
in equilibrium, the disorder potential values must lie be- 
tween ±{l/2)C~l^. The magnitude of this quantity equals 
i?**' in the onsite limit and decreases as the strength 
and range of the Coulomb interactions increases, see Ap- 
pendix II. Additionally when the array is in equilibrium 
state, the energy to hop between all pairs of adjacent 
sites must be greater than zero. From the disorder 
potential differences are restricted between ±E'^^^. 
El^^ equals 2i?**' in the onsite limit and decreases with 
decreasing d/r'^^\ as shown in Fig. 8. 

Fig. 1(a) and (b) compare P{(j)f') and P{^f') for 
arrays with purely onsite interactions {C^_lj — 0) with 
arrays with long range interactions (C^J'j ^ 0) at two 
spacings, d/r^'^^ = 0.5 and d/r"^"^ = 10. In all cases, 
and are calculated for arrays with 50 islands in be- 
tween two grounded leads. The histograms average the 
values of the potentials of all islands and the values of the 
potential drops between all adjacent islands over many 
realizations of disorder (0(> 10^)). The smaller spacing, 
^jj-isi _ Q 5^ jg typical of chemically assembled nanopar- 
ticle arrays. The larger spacing, d/r"^^^ — 10, is atypical 
of arrays in most experiments but is mainly included as 
a pedagogical example because it has interactions among 
islands that are finite yet comparable to the onsite case 
that is often used to describe experiments^i^i^. Ar- 
rays recently synthesized^ have large rf/r*"' values. In 
Fig. 1(a), P(0f'*) is a constant between ±i5J*' for the on- 
site case. As the range of interactions increases (decreas- 
ing d/r*"'), the width of P(0f'*) decreases because the 
disorder potential values are bound by ±0.5(7^^^. Increas- 
ing Coulomb interactions also increases (decreases) the 
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FIG. 1; Probability distributions of disorder potentials (j> " in 
(a), and disorder potential diflterences in (b), due to dis- 
order for 50-island arrays with purely onsite interactions and 
with long range interactions at two spacings: d/r'"' = 0.5 
and 10. The onsite distribution is described by Eq.JTS]) 
with A^i^Ax = £^'=7^;°""*'= = where E'/' = 1/(2C"')- 

For all the cases, P((^*'') (P($'*''')) is finite valued for 0'*'" 
(<1>'^") between ±0.5C-^ {±E^-^) which decreases with de- 
creasing spacing. As the spacing decreases, the probability of 
having and values close to zero increases. Vertical 
lines are included as guidelines to emphasize the edges of the 
distributions. 



probability of small (large) values of |(/)**| . In Fig. 1(b), 
the onsite P($f*'*) distribution is given by p^ . Simi- 
lar to the trends in Fig. 1(a), as the range of Coulomb 
interactions increases, the width of the distribution de- 
creases and the probability of small (large) |$**| val- 
ues increases (decreases). The increased probabilities of 
small l^f'"! are due to Coulomb correlations that make 
it more likely for the disorder potentials of neighboring 
islands to have similar values. See Fig. 2. Increasing 
the range of Coulomb interactions leads to a greater rel- 
ative reduction in the width of P{^f'') than Picj^f") be- 
cause the former are bound by ^E^^^ whereas the latter 
are bound by ±0.5C,^^. In the onsite case, E'i~^ equals 
2E]!'^ for all junctions between two islands and increasing 
Coulomb interactions can reduce Ef~^ significantly due 
to an decrease (increase) in C~l^ (C'~jl|-i). See dH). 

Our results for P{^f^'') differ to some extent from those 
by Elteto et al^, calculated with an inverse capacitance 
matrix C^^ that is finite only for nearest neighbors and 
charge disorder modelled by a set of stationary quenched 
charges. Elteto et al^ distributions are bound by ±Cj7^ 
instead of by E^~^ due to the lack of correlations in their 
quenched disorder model. We permit the interactions 
among the screening charges to determine whether or 
not the disorder potentials are correlated. 




Site Index (k) 



FIG. 2: Comparison of < ^i't^-f " > normalized by l^i'p for 
50-island arrays with onsite (top plot) versus long range (bot- 
tom plot) Coulomb interactions. In the presence of long-range 
interactions between the charges, the values of the disordered 
potentials are also correlated. CJ^25 normalized by C2525 ^ 
included in the long range case to show that correlations in 
the disorder potentials are related to, but decay faster than 
the elements. 

In Fig. 2, we plot < (f)'^'^ > to show how interac- 
tions among charges affect the correlations among disor- 
der potentials. In the onsite case, correlations are fi- 
nite only if i = j. In this case, the disorder poten- 
tials of different islands are uncorrelated. In the case 
of long-range interactions with d/r*^' = 0.5, correlations 
are maximal when i — j, but they do not vanish for 
j. < 0f > is finite for at least |i - fc| < 3 - 4. 
The correlation of disorder decays faster than the inter- 
actions as shown in the figure. The correlations between 
and its nearest neighbors make it more likely 
for the disorder potential differences to have small 
magnitudes. 



IV. THRESHOLD 

The threshold voltage is the minimum voltage which 
allows the flow of current. In this section we analyze the 
threshold voltage of one-dimensional arrays with long- 
range interactions for both clean and disordered systems. 
As in the short-range case the threshold is completely 
determined by energetic considerations and is indepen- 
dent of junction resistance disorder. In the clean case 
we flnd that the threshold equals the minimum voltage 
necessary to create an electron-hole pair and increases 
with the number of particles at a rate which depends on 
d/r*^'. For charge disordered arrays the average thresh- 
old is reduced compared to the onsite interactions case. 
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and increases linearly with the number of islands, with 
a slope that decreases with decreasing d/r''*', i.e. with 
increasing the range of the interaction. 

When the interaction strength between charges at the 
nanoparticle and those at the electrodes does not vanish 
for any particle the polarization potential drop at every 
junction is finite. In a clean array, the potential gradient 
created by this polarization potential drop allows current 
once an electron-hole pair is created, opposite to what 
happens in the onsite case. As a result, the threshold 
voltage equals the mimimum voltage which allows the 
creation of an electron-hole pair. 

The cost in energy to create an electron-hole pair 
in junction i in an uncharged array (i.e. an array 
in which the nanoparticles have no excess charges) is 



AE = Ef 



-h 



AfV. We can define a junction depen- 



dent threshold voltage for creating an electron-hole pair 
V^"'" = E^-'^/Af. In the onsite limit V,^"" is finite 
only at one or both contact junctions and infinite at the 
bulk, but with long-range interactions V^"^^'" is finite at 
every junction. For those cases analyzed, we have found 
that due to the smaller value of the excitonic energy and 
the larger potential drop Vj^^ is smallest at the contact 
junctions and the threshold voltage is controlled by them. 

Figs. 3(a) and (b) show the dependence of the thresh- 
old voltage Vt of clean, symmetrically biased arrays with 
long-range interactions on the number of islands in the 
array N and on the spacing between array sites, d/r'*'. 
The threshold voltage is determined by those factors that 
define the polarization potential drops across the con- 
tacts, A" and A'^^^. Vt increases with increasing N 
because the fraction of the polarization potential which 
drops across the contact junctions decreases as N in- 
creases. As the spacing between the leads increases, the 
polarization potential drop across each contact decreases 
until eventually it reaches a minimum value at which the 
polarization of each contact is only due to the interaction 
of each contact with the lead closest to it. As a result, Vt 
increases sublinearly with increasing N and eventually 
saturates. For N and d/r"^^^ large enough that the polar- 
ization potential drop across the contacts is not strongly 
influenced by interactions with the opposite lead, de- 
creasing the array spacing decreases the polarization po- 
tential drop across the contact junctions and the thresh- 
old increases. For N and d/r^"^ small enough that both 
leads strongly influence the polarization of both contact 
junctions, decreasing the spacing increases the polariza- 
tion potential drop across the leads and the threshold 
decreases. The potential drops and threshold can be es- 
timated by using an unscreened model for the inverse 
capacitance elements associated with the leads, C^^ and 

Ci"7V-(-i- These estimates are included as dashed lines in 
Figs. 3(a) and (b). 

As shown in the inset of Fig. 3(b) Vt changes smoothly 
with a, contrary to the peak-valley structure found in the 
onsite interaction case^. The threshold voltage depends 
less on a as iV increases and as d/r'*^' decreases because 



30 



20 



> 



10 







d/r" 


=10 


d/r" 


'=1 


d/r" 


=0.5 



9/' 



_L 



_L 



_L 



1 

r 0.75 Y- / - 



0.5 



_L 



V N=50 



0.5 1 
a 



N=5 



50 



40 



30 



20 40 60 2 4 6 8 10 



A 

>^ 20 
V 

10 



10 




^ ' 

^ " 


to 




^ ■ 

^ ■ 


1 


1 1 III 


^ " 


10 


^ ■ 

^' 


1 


' 1 , 1 


1,1, 



10 



20 30 
N 



40 



50 



FIG. 3: (a) and main figure in (b); With symbols it is respec- 
tively plotted the threshold voltage Vt of symmetrically bi- 
ased arrays a = 0.5, with no disorder as a function of number 
of islands N and of array spacing d/r"'' for the inverse ca- 
pacitances calculated as described in appendix I. As decribed 
in the text, the threshold voltage of clean arrays is controlled 
by the value Af at the contact junctions. The dashed lines 
are estimates of the threshold voltage that use a inter- 
action to approximate the polarization potential drops across 
the contact junctions A" and A^_|_i. Inset in (b): Thresh- 
old voltage of clean arrays for several array parameters as 
a function of a normalized to the value for a symmetrically 
biased array. From top to bottom d/r'"'' = 0.5, A'^ = 50; 



0.5, A = 20; and d/r'- 



10, A = 50. (c) Main 



figure: Average threshold voltage of disordered arrays versus 
the array length at three different array spacings. The solid 
line shows the dependence of the average threshold on array 
length in the limit of onsite interactions. The inset shows the 
root mean square deviation of the threshold voltage distribu- 
tion. For small d/r'"' it deviates from the N^^^ dependence 
found in the onsite case (in solid line). In inset and main fig- 
ure dashed-dot lines are a guide to the eye and same legend 
as in (a) applies 



the applied voltage drops more homogenously across the 
array junctions. This dependence disappears completely 
if the polarization potential drops linearly. In this last 
case the threshold voltage of clean arrays would be Vt = 
(A^ -I- 1 )_£**'. This threshold is linear in the number of 
particles, as in the onsite case, but its origin of linearity 
and its slope differs in both cases. 
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In the charge disordered onsite case every up-step in 
the disorder potential prevents current flow and has to 
be compensated by a charge gradient. As a result, the 
threshold voltage is controlled by the number of up-steps, 
but this is not the case for long-range interactions. In the 
long-range case the situation is more complex. Due to the 
finite polarization potential drop at the inner junctions 
the number of junctions which prevent the flow of current 
is reduced, compared to the onsite case. Up-steps in the 
disorder potential drop can localize the charge only if its 
value is larger than the polarization potential drop at 
the same junction. In some cases in which there is finite 
polarization drop at a junction, that is slightly smaller 
than the energy cost for tunneling it is possible that a 
small increase in the voltage in the electrodes permits 
the tunneling. Increasing the voltage above the minimum 
bias voltage which permits the generation of electron- 
hole pairs will, in some cases, but not always, result in a 
negative potential drop at the given junction allowing the 
flow of charges. But, quite often, the entrance of more 
charges and the creation of a charging potential gradient 
will be required in a similar way as it happened in the 
onsite case. Note that due to the polarization voltage 
drop at the bulk, the threshold junction can be other 
than the contact ones. Moreover, the interaction between 
charges in different islands decreases the energy for the 
entrance of a charges with opposite sign and increases 
the one for the entrance of charges with the same sign. 
This effect was attributed to an attraction (repulsion) 
between the injected soliton and an antisoliton (soliton) 
on the array by Likharev and coauthors^. Accumulation 
of charges in the array increases the threshold voltage. 
On the other hand a value of at the contact islands 
favourable for the entrance of charge unto the array can 
decrease it, as the polarization potential drop to allow 
entrance of charge is smaller. Both mechanisms compete 
to determine Vr- For large d/r^^^ the accumulation of 
charges is more important as the voltage drop at the bulk 
junctions is small and on average Vt of large arrays will 
increase compared to the clean array threshold voltage. 
On the contrary the second effect can be more important 
for small d/r^^^ . 



Numerically we have found a linear dependence of the 
threshold voltage on the number of particles in the ar- 
ray, see Fig. 3(c). Decreasing the array spacing decreases 
the average thresholds below the threshold values of the 
arrays in the onsite limit. The expected behavior of the 
average threshold value as compared to the clean limit, 
discussed above, is found. Only at the largest array 
spacmg (d/r^^' = 10) studied do we recover the depen- 
dence of the fluctuations of the threshold voltage on array 
length predicted by Middleton and Wingreen^i, see inset 
of Fig. 3(c). 



V. FLOW OF CURRENT 

In this section we discuss the three main voltage 
regimes which can be distinguished in the I-V curve. At 
voltages very close to the threshold, we show that the cur- 
rent is linear in {V — Vt)- Contrary to the short-range 
interaction case, the slope of this linearity depends, not 
only on the junction resistance and on a, but also on the 
number of islands and on d/r'*'. For given a, N , d/r'^^^ 
and {Ri], in the presence of charge disorder, the slope of 
the I-V curve close to threshold can depend on the po- 
tential disorder conflguration. At intermediate voltages 
steps in the current are smoothed compared to the on- 
site interaction case. A linear I-V with an offset voltage, 
closely related to the one found for short range interac- 
tions, characterizes the high voltage regime. 

A. Voltages close to threshold 

For the onsite interaction case we recently resolved the 
controversy^^i^i^i^ on the power law dependence of cur- 
rent on iy — Vt), showing that very close to threshold it 
is linear--. Such a linearity is due to the bottle-neck char- 
acter of one of the junctions and the linear dependence of 
the energy for tunneling on the bias voltage. These two 
assumptions remain valid for long-range interactions, so 
the linear {V — Vt) dependence is also found in this case. 
In the case of clean arrays, the threshold and low-voltage 
bottle-neck for current are found at the contact junctions 
and 

I^^§^{V-Vt) (15) 

Except for a — 1/2, for which Aj^'^/i?! and h^Jj"]^^/ Rn+i 
have to be added in the expression for the slope. Depen- 
dence in TV, d/r''^^ and a via the dependence of A"^_|_]^ 
on these parameters is found as seen in Fig. 4 (a) and 
(b). The value of Af which appears in the expression 
of Vt is the same that controls the linearity of the IV 
curves very close to theshold. The behavior of the slope 
of / vs (y — Vt) with a, N and d/r^^^ is opposite to 
the one of Vt- For a different from 1/2, increasing N 
and decreasing d/r'''' decreases the slope because these 
changes reduce the fraction of the polarization potential 
that drops across the contact junctions. Biasing the array 
in a more assymetric way (increasing ( | a — 1/2 1 ) changes 
the slope as the fraction of the polarization potential that 
drops across the contact junction that serves as a bottle- 
neck at small voltages is modified. The slope does also 
depend on the junction resistances, similar to the depen- 
dence found for onsite interactions^^, see Fig. 4(c). 

The charge disordered case is to some extent different. 
The above threshold bottle-neck (and below threshold 
current-blocking) junction is not necessarily either of the 
contact junctions, so the slope of the linear dependence 
can be controlled by Af with i ^ 1, -I- 1. The junction 
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rectly modified by changes in the bias potential only at 
the contact junctions, at least for a = 1/2. The slope is 
generally larger (smaller) when the bottleneck junction 
lies closer to the edge (middle) of the array. Changing a 
also modifies the slope of the disordered case, not shown. 
This modification can be due to the change of Af with 
or without a change in the bottleneck junction. 



B. Intermediate voltage regime 

As in the onsite case^S the linearity of the current 
disappears when the bottle-neck description stops being 
valid. This happens at very small values of {V — Vr) ^ 
lO^^EIfK It is easy to show how this situation leads to 
sublinear behavior. Assume that the transport happens 
through a sequence of iV-t- 1 tunneling processes and con- 
sider a bottle-neck process with rate Ti — R^^AfV with 
V = V — Vt and another process in the sequence with 
rate Tj — R~^{Ej + A'jV). Here Ej is the gain in energy 
of the second process a.t V = Vr- If these two processes 
have rates much smaller than the rest of processes in the 
sequence, the current can be approximated by 



(V-VJ/E^ 



FIG. 4: I-V curves of ordered arrays at voltages very close to 
threshold. The insets show the derivatives of the I-V curves. 
Similar to what was found for onsite interactions very close to 
threshold the I-V is linear. The linear regime ends at voltages 
Y — Yj, ~ 10~^El"'' , much lower than the values used in exper- 
iments to check the power-law dependence close to threshold, 
(a) shows how varying the length and spacing of symmetri- 
cally biased arrays alfects the slopes of the linear regimes, 
due to the change in the polarization potential drop factors 
Af and correspondingly the fraction of potential which drops 
in the junction which acts as bottle-neck, (b) In a similar way, 
the slope depends on how symmetrically is the array biased. 
(c)plots the I-V curves and derivatives corresponding to clean 
arrays with equal length, d/r"' and a but different junction 
resistances. Resistances in these plots vary randomly in be- 
tween (5 - 11)7?T, (8 - 21) Rt and 23 - 83Rt in top, middle 
and lower curves. (d)I-V curves corresponding to three differ- 
ent realizations of charge-disordered arrays with all junction 
resistances equal TV = 50 and a = 0.5. Contrary to what was 
found for onsite interactions the slope of the low-voltage lin- 
ear current can be dilferent for different arrays with the same 
nominal parameters if there is charge disorder. This reflects 
that the bottle-neck is not necessarily a contact junction but 
can be a junction in the bulk. 



that acts as the bottle-neck depends on the particular 
disorder configuration. This is shown in Fig. 4(d), which 
plots the current and its derivative with respect to volt- 
age for several configurations of disorder corresponding 
to the same value of d/r**', N and a and the same junc- 
tion resistances. A similar dependence of the slope on 
the charge disorder configuration was not present in the 
short-range case, as the energy gain for tunneling is di- 



i?r^A?y 



R-^AfV 1 



r.-\'ejIajv) 
R7^A?v\ 



R-'E. 



(16) 



The slope of the current and the lost of linearity de- 
pends on the resistance of the junctions, as was also seen 
in the onsite case^. In the clean long-range case com- 
paring the values of Af the linear behavior lasts longer in 
shorter arrays, smaller d/r'*' and smaller a, as in these 
cases the values of A" are more homogenous through- 
out the array. The disordered long-range case is more 
complex. Due to the non-homogeneous increase in polar- 
ization voltage drop a junction which has a small energy 
gain can increase this gain more than other junctions 
when the applied voltage increases and the dependence 
of the slope with the array parameters is not so easily 
predicted. 

If the value of the voltage is increased further several 
tunneling processes are energetically allowed at each step 
in a sequence and the discussion of transport becomes 
more complex due to the multiple choices available and 
the polarization potential drop at the bulk junctions. 
Above the linear regime there is a region of smoothed 
steps in the I-V curve. Decreasing d/r*'*' smooths the 
steps and for small d/r*'*' they are hardly distinguish- 
able. This behavior is seen in Fig. 5(a) which compares 
the onset of current, at voltages not extremely close to 
threshold, for several array parameters. For clarity the 
curves have been plotted as a function of {V — Vt)- The 
staircase profile differs in all these cases. The top curve 
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FIG. 5: In (a) and (b) the I-V characteristics show the 
Coulomb staircase and correspond to clean arrays with ho- 
mogeneous resistances, (a) shows five curves for diflerent in- 
teraction ranges and lengths. From top to bottom, the first 
four curves correspond to a = 0.5 and respectively to A'' = 50 
with onsite interaction; = 50 with d/r'°' = 0.5; = 30 
with d/r'"' = 10; iV = 50 with d/r"' = 10. The lowest curve 
corresponds to a A'' = 50 and a = 1 and d/r"' = 10. The 
threshold voltage has been substracted. Vr equals 98, 43, 14, 
16 and 9.7 respectively. In the onsite case it is necessary to 
create a charge gradient at every junction to allow current, 
but once charge can reach the opposite electrode it fiows eas- 
ily and the TV shows a large jump close to threshold. For 
d/r"' = 10 a charge gradient is not created and bulk junc- 
tions slow the current flow producing the staircase structure. 
The step width is smaller for forward bias a = 1 than for 
symmetric bias, but contrary to the behavior found in the 
Coulomb staircase for onsite interactions, with long-range in- 
teractions the step width is not flxed. The steps are washed 
out for d/r"' = 0.5 due to a more homogeneous polarization 
potential drop through the array, (b) I-V curves for disor- 
dered arrays with d/r"' = 0.5, 10. The clean case is included 
for comparison, (c) I-V in a large voltage scale for clean ar- 
rays with homogeneous resistances , from top to bottom solid 
lines correspond to d/r"' = 0.5, d/r"' = 10 and onsite inter- 
actions. The dashed lines give the asymptotic predictions for 
d/r"' = 0.5 and 10. At high voltages all the curves have the 
same slope given by R^um- The offset voltages differ as the 
excitonic energies do. In spite of the very different thresh- 
old voltage and the different dependence on voltage close to 
threshold, the d/r"' = 10 I-V curve differs little from the 
onsite one at high voltages, (d) TV curves in a large voltage 
scale corresponding to d/r*"' = 0.5 symmetrically biased ar- 
rays. From top to bottom A'^ = 30 and A' = 50 disordered 
arrays with homogeneous resistances and a A^ = 50 clean ar- 
ray with the first resistance ten times larger than the other 
ones. Dashed lines give the asymptotic predictions. The slope 
of both A' = 50 curves differ, but their offset voltages are the 



corresponds to an = 50 array with short-range inter- 
actions, previously studied^^. In this case to allow cur- 
rent flow a charge gradient at each bulk junction has to 
be created, but once charge can enter the array it flows 
easily through it. This is reflected in the sharp onset 
of the current close to threshold. The steps at higher 
voltages are just barely visible at this scale. Once the 
polarization voltage at each junction is finite it is not 
necessary to create a charge gradient at the inner junc- 
tions and the steps' shape is modified. The three bot- 
tom curves correspond to d/r*'*' = 10. Several features 
can be appreciated in these curves. The step shape is 
clearly seen. As the potential drop at the inner junc- 
tion is small, bulk junctions control the flow of currents 
at each plateau and their slope is small. The slope is 
larger for shorter arrays, and the step width depends on 
the value of a. The dependence of the step width on a 
was also found for onsite interactions, but contrary to the 
short range case, for finite d/r'^'^^ the step width is not a 
constant throughout the curve as the presence of charges 
in the array influences the energy cost to add charges 
from the electrodes, to the first or last island. As also 
seen in this figure, for small d/r"*' the polarization po- 
tential drop at the inner junction is larger and similar to 
what happened in the onsite case (but for reasons to some 
extent different) once the charge enters the array it can 
flow easily. In (b) we can see the different I-V curves in 
clean and disordered arrays. Specially interesting is the 
disordered d/r*'*' = 0.5 TV characteristic. It looks su- 
perlinear, similar to what would be found if a power-law 
larger than unity is present at these voltages. We have 
observed that this approximate superlinear type depen- 
dence is common in disordered arrays with small d/r**'. 
If experimentally the power-law behavior expected close 
to threshold is measured at these voltages (larger than 
those at which the linear behavior is predicted) the ex- 
ponent of the power-law could be erroneously assigned a 
value larger than one. 



C. Linear behavior at high voltages 

At very high voltages, in the onsite interaction case we 
showed^^ that the current can be approximated by the 
asymptotic expression 



1 



N+l 



^asympt 



with Rsum = X^ila^ ^i- The arguments, based on a uni- 
form tunneling rate through all the junctions in this volt- 
age regime, which led to this expression remain valid in 
the long-range case, with the only change of the quanti- 
tative value of the excitonic energies E^~'^- The slope of 
the current does not depend on the range of the inter- 
action or the presence of charge disorder but the offset 
voltage at which the asymptotic expression cuts the zero 
current axis, does^i^i^. Confirmation of the validity 
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FIG. 6: Average voltage drop close to threshold for A'^ = 50 
arrays with different parameters. Both (a) and (b) correspond 
to clean arrays. Main figures are for a = 0.5 and long-range 
interaction. The average potential drop essentially equals the 
polarization potential for each value of d/r'"', which is plot- 
ted as filled small dots in (b) for comparison. This behavior 
contrasts with the potential due to the charge gradient which 
has to be created in the onsite case, shown in the inset of (a). 
As shown by comparing main figure in (a) with the inset in 
(b) a change in the value of a modifies Ai and correspondingly 
the potential drop through the array. 



of pT)l is seen in Figs. 5(c) and (d), where numerical 
results are compared with the asymptotic behavior pre- 
dicted by pT|) . In Fig. 5(c) the I-Vs for symmetrically 
biased arrays with iV = 50 nanoparticles and different in- 
teraction range are plotted. At high voltages slopes are 
equal but the offset voltage to which they extrapolate it 
is not. Note the difference between the value of threshold 
and the one of the offset. In particular the d/r*"' = 10 
curve has a smaller threshold and larger offset than the 
^j^^si _ Q pjg^ g^j-j influence of the number of 
junctions and their resistance in the high- voltage current 
is reported and shown to be in good agreement with the 
approximate prediction. 



VI. VOLTAGE DROP 

In this section we analyze the effect of long-range in- 
teractions in the voltage drop through the array. We 
differentiate the same three regimes as in previous sec- 
tion. Differences found with respect to the onsite case, 
previously analyzed^^ are mainly due to the polarization 
potential drop at the bulk junctions, which vanishes in 
the onsite limit but is finite when interactions between 
the charges at the islands and those at the electrodes are 
long-range. 

As discussed in section IV, in a clean array, Vr is given 
by the minimum bias voltage which allows the creation 
of an electron-hole pair. Contrary to the onsite case^ 
there is no charge accumulation in the array. It is thus 
expected that very close to Vr the voltage drop reflects 
the polarization drop A"y at each junction. The po- 
tential drop distribution thus depend on the length of 



the array iV, the bias asymmetry a and on the range 
of the interactions djr'^^^. This dependence is confirmed 
in Figs. 6(a) and (b) where the potential drop is plot- 
ted for different array parameters. The value of A" is 
included for comparison in Fig. 6(b). The voltage drop 
is very different from the one found to the onsite case 
(included in the inset of Fig. 6(a)), where in the bulk it 
is due to charge accumulation at the islands. The depen- 
dence on the value of the resistance is extremely weak 
even once the polarization potential drop is substracted 
(not shown) and not observable, except if the difference 
in the value of resistances is very large. In the disordered 
case with long-range interaction it is possible that once 
charge is allowed to enter the array it can flow. Then the 
average potential drop is approximately the sum of the 
disorder potential and the polarization potential. In gen- 
eral, when this happens the threshold voltage is smaller 
than the one in the clean case as the disorder potential 
reduces the polarization potential drop necessary at at 
least one of the contact junctions. But for large d/r*'*' is 
more probable that one or more charges remain stacked 
in the array, similarly to the case of onsite interactions 
and the charge potential due to these stacked charges has 
to be added. 

At intermediate voltages, in the Coulomb staircase 
regime, we saw in the onsite case^ that the voltage drop 
through the array shows an oscillatory behavior, with 
the number of maxima increasing from step to step. A 
similar behavior is found in Fig. 7(a) corresponding to 
a clean array and d/r''*' = 10. In Fig. 7(b) for all the 
voltages plotted the number of maxima is two, and their 
amplitude decreases until at the largest voltage oscilla- 
tions cannot be discerned. Comparing with Fig. 5 one 
realizes that the step number has not changed. The I-V 
curves reaches the high-voltage regime without showing 
stepwise behavior. 

At high voltages, the potential drop is qualitatively 
similar to the one found in the onsite case^. The volt- 
age drops linearly only after subtracting from each junc- 
tion the excitonic energy. The sum of the excitonic ener- 
gies results in the offset voltage. The current is equal to 
iy — Vof fset)/ Rsum and the corresponding voltage drop 
at each junction is 



E' 



(18) 



which satisfies V = = +J2iEr'' = 

I Rsum + Vof fset- Eq. p8|) is valid for both ordered and 
disordered arrays. Some examples of this behavior are 
shown in are shown in Fig. 7. In Fig. 7(c) we can see that 
as expected, in the absence of resistance disorder, once 
the excitonic energy is substracted the potential drops 
homogeneously through the array even if there is charge 
disorder, while it is proportional to the resistance value 
when resistances are not homogeneous, as in Fig. 7(d). 
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FIG. 7: Average potential drop at the array junctions at in- 
termediate (upper figures) and high voltages (lower figures) 
for Q = 0.5 and N = 50. Upper figures correspond to clean 
arrays. Curves in (a) are for d/r"'' = 10 and (from top to 
bottom) V = 20,32,44, 56-E^'''. The oscillations in average 
potential found in the onsite case at intermediate bias volt- 
ages are still present, and the number of maxima increases 
in two when going to a higher step. In (b) d/r"' — 0.5 and 
V = 46, 60, 70, 80, QOEl"' . The number of maxima has not in- 
creased in this range of voltages but oscillations are smoothed 
with increasing bia voltage. In spite of the large value of the 
voltage the current is still in the first step of the Coulomb 
staircase, see Fig. 5(a). The potential drop at V = 90-E"' re- 
sembles the one at high voltages, homogeneous except by the 
excitonic term responsible of the offset. In (c) and (d) the ex- 
citonic energy has been substracted from the average potential 
drop at high voltages. Once this term is substracted the av- 
erage potential drop is completely homogeneous through the 
array in (c) where there is charge disorder, and all resistances 
are equal but not in (d) which corresponds to a clean array 
but with resistances which vary randomly between (5 — 11)-Rt 



VII. SUMMARY 

In this paper we have analyzed the effect of the long- 
range interaction in the transport through ordered and 
disordered nanoparticle arrays. To this end we have in- 
troduced a model which allows a discussion of the rele- 
vant quantities in determining these transport properties. 
We have introduced the concept of polarization potential 
drop, which results from the interaction between charges 
at the electrodes and at the islands. In the model used 
we take into account that electrodes are not ideal voltage 
sources. Their potential fluctuates in response to tunnel- 



ing processes (but its effect in the numerical results is 
very small). 

We have studied how the proximity of the nanopar- 
ticles modifies the 1/r-interaction due to screening. To 
determine this screening we have developed two meth- 
ods which allow us to calculate the inverse capacitance 
matrix of the system under study, see Appendix I. As 
discussed in Appendix II, and shown in Fig. 8(a), for the 
case of metallic nanoparticles, here considered, the effect 
of screening starts to be relevant when d/r^^^ ~ 1 — 2, 
there is no divergency in the value of C~j^^ at small djr^^^ 
values, but capacitance values and their inverses saturate 
at finite values. As discussed previously by Matsuoka 
and Likharev^^ for cylindrical nanoparticles, the interac- 
tion between charges is reduced only when the nanopar- 
ticles are very close, at larger distances the interaction 
increases and approaches the 1/r law from above, result- 
ing in a bump in the interaction potential as compared 
to Coulomb law, see Fig. 8(b). We have related this 
anti-screening effect with the dipolar charges induced in 
the conductors. 

Long-range interactions screen the disorder potential 
and induce correlations between the values of (i** at 



different islands. These correlations are related to, but 
smaller than the inverse capacitance matrix elements. 
The distribution of the island and junction potentials is 
modified in comparison to the one found for onsite in- 
teractions. These effects are discussed in section III, and 
are shown in Figs. 1 and 2. 

As discussed in section IV, as in the onsite case, the 
current is blocked up to a threshold value Vt- This 
threshold value is independent of the resistances of the 
junctions. In clean arrays, Vr is the minimum volt- 
age that allows the creation of an electron-hole pair, see 
Figs. 3(a) and (b). This behavior is opposite to the one 
found in the onsite case where a charge gradient at the 
junctions has to be created to allow the flow of current. 
With long-range interactions, the polarization potential 
drop across the array creates a potential gradient which 
facilitates charge flow. In the disordered case, two ef- 
fects compete that can increase or decrease the threshold 
voltage compared to the clean case. Charge accumula- 
tion can be induced by up-steps in the disorder potential, 
increasing Vr, and the disorder potential distribution can 
reduce the energy to create an electron-hole pair decreas- 
ing it. The latest effect dominates for small d/r*'*', see 
Fig. 3(c). 

The current varies linearly with respect voltage close 
to threshold, in the bottle-neck regime. This is discussed 
in section V and plotted in Fig. 4. The slope depends, 
as the threshold voltage does on the polarization poten- 
tial drop, and also on the resistance of the bottle-neck 
junction. At intermediate voltages, we find steps in the 
current, but these are smoothed as compared to the on- 
site case, see Fig. 5(a) and (b). Contrary to the onsite 
case, due to the interaction between the charges in dif- 
ferent islands the width of the steps is not fixed. A linear 
dependence of current on voltage is recovered for large 
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biases, but as in the onsite case, current extrapolates to 
zero at a finite offset voltage (Fig. 5(c) and (d)). The off- 
set voltage is given by the sum of the excitonic energies 
of all the junctions, and its value depends on the range 
of the interaction, but as in the case of onsite interac- 
tions the slope of the current is given by the sum of the 
resistances in series and is independent of d/r"^^^ . 

The voltage drop through the array close to thresh- 
old, reflects the polarization contribution A" due to the 
electrodes, as shown in Fig. 6. In the onsite case, the 
potential drop in this bottle-neck regime was due to the 
charge accumulation at the islands, necessary to create a 
potential gradient through the array. 

Shown in Fig. 7, at intermediate voltages an oscillatory 
voltage drop through the array, similar to the one found 
for short range interactions is found for large values of 
d/r*'*'. For small d/r*'*' there is some remnant of this 
behavior but the amplitude of the oscillations can vanish 
while still being in the first step of the Coulomb staircase, 
which is not well defined in this case. At large voltages, 
the potential drop is analogous to the one found in the 
onsite case. 

In conclusion, we have found that the long-range char- 
acter of the interaction modifies the I-V characteristics of 
one-dimensional arrays for both ordered and disordered 
arrays. This modification is greatly due to the finite po- 
larization potential drop at the bulk junctions, due to 
the interactions between charges at the islands and at 
the electrodes. Differences in transport as a function of 
the range of the interaction are larger at low voltages. 
In both the long-range and onsite cases the current is 
blocked up to a threshold voltage Vt and the current de- 
pends linearly on for very small {V — Vt), but the value 
of both Vt and the slope of the current are different in 
both cases, as the mechanism to create a potential gra- 
dient through the array differs. At high voltages, for 
given values of the excitonic energies the transport in 
both cases is analogous. 



VIII. APPENDIX I: METHODS TO COMPUTE 
THE CAPACITANCE MATRIX 

In this appendix we discuss the two methods that we 
have used to calculate the interaction strength C~^. In 
the first one we determine the capacitance matrix using 
the method of images. This is an iterative method which 
a priori can be used to determine the capacitance ma- 
trix Cq/3 for any geometric configuration of spheres, so 
it is valid also for two-dimensional arrays, for example. 
Although the algorithm for generating images is straight- 
forward, the number of images required to calculate Caf3 
makes the numerical implementation of this technique 
nontrivial. While computer memory problems can be 
solved, the computation time is too large to tackle those 
cases with very large arrays and electrodes and small dis- 
tance between conductors. On the other hand, the ca- 
pacitance matrix is calculated to a given accuracy. Small 



errors in Cap can be enhanced and uncontrolled in C^i^ . 
In the second method the interaction matrix C~l is cal- 
culated directly taking into account the symmetry of the 
system and the properties of spherical harmonics. It re- 
quires the inversion of a matrix which can be quite large. 
It is specially useful and fast for systems with azimuthal 
symmetry as the one considered here. Results obtained 
with both methods are in extremely good agreement. 



A. Image Charges Method 

The method of images is based on the relation between 
charges and potentials in capacitively coupled conduc- 
tors. The charge Qa induced on a conductor in the pres- 
ence of K equipotentials at potentials Vg is given by the 
capacitance matrix Cap 

K 

/3=1 

The inverse capacitance matrix which enters the free en- 
ergy IT]) is the inverse of the capacitance matrix Cap ■ We 
have calculated the capacitance matrix using some prop- 
erties of spheres and the fact that the charge in a con- 
ductor a produced by a unit potential in a conductor (3 is 
equal to the capacitance matrix element Cap- We have 
obtained C~p inverting Cap- The method of images^ 
is the placement of imaginary charges inside the spheres 
at positions that make the potential everywhere on the 
surface of the conductor equal a constant. 

To determine the positions of the image charges, we 
exploit two properties of spheres. First, the center of the 
sphere is equidistant from all points on the surface of the 
sphere. Using this property, the surface of an sphere of 
radius R can be set to a potential V by placing an image 
at the center of the sphere of charge q = VR. Second, for 
every point outside a sphere there is a point inside the 
sphere for which the ratio of the distances between these 
two points and any point on the surface of the sphere 
is a constant. From here it follows that if a charge 
is located at the outside point, at a distance dc from the 
center, an image charge qi placed at the inside a distance 
R? jdc from the center in the radial line, with charge 

qi - -qRj (20) 

will set the potential to zero everywhere on the surface of 
the sphere. We determine the {N -I- 2) x (iV + 2) capaci- 
tance matrix, column by column, by determining the set 
of image charges that sets the potentials of the spheres 
to Va — <5a/3- The capacitance matrix elements Cap are 
given by the sum of all the charges in sphere a. To set 
the potential of the /3 sphere, with radius Rp to one, 
we place a charge with magnitude Rp at the center of 
this sphere xp. The remaining spheres are grounded by 
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placing images inside each sphere with charges 

Ui,q 



(21) 



The multipolar charge induced is the one which mini- 
mizes the energy. Separating the monopolar contribu- 
tion {I, m — 0) in the expression of the free energy, and 
minimizing the latter with respect to Qf„i, we obtain 



at positions 



R„ 



(22) 



Here q°^'^ and Xg^,^ are the value and the position of the 
charge which creates the inhomogeneous potential that 
we want to compensate and Ru and the radius 

and position of the center of the sphere to which we add 
the image charge q,y. These image charges are added to 
all the spheres except the one in which qoid is placed. 
The charges that have been added generate new inhomo- 
geneous potentials at the rest of the spheres and have to 
be compensated following the same method. This pro- 
cess repeats iteratively for all the charges added to all 
the spheres. During each iteration n, the number of new 
images required to compensate the potential of the other 
spheres approximately equals {N + l)". We eliminate 
some of the images by discarding images with a magni- 
tude that is smaller than a suitable cutoff value, qcutoff- 
We required qcutof / to be small enough that the rela- 
tive diffferences between the matrix elements generated 
with the cutoff value qcutof f and by a larger cutoff value 



^cutoff 



10* 



cutoff 



are less than one percent. 



B. High-order multipoles method 

Following Wehrli et a*'', the energy of the system, 
given by ([T]) can be rewritten in terms of the higher- 
order multipolar charges induced by the charges on the 
conductors as 



2 



a,* ^0.(3 
Lm LmA' 



(23) 



Here Greek indices denote the conductors, / and I' de- 
note the order of the multipole, and m — —I, . . . ,1 and 
m' = — Z', . . . ,1 denote the azimuthal number. This ma- 
trix G is hermitian with respect to the exchange oia,l,m 
and P,l',m'. Using the linear response form for the in- 
duced multipoles, the higher-order multipolar charges, 
QTm^ can be expressed in terms of the (monopolar) 
charges on the conductors — Q2q, as 



Qa 



(26) 



Here A — l,m and / 7^ 0, correspondingly B, and the 
equation is written in vectorial and matrix notation. 
Once this expression is substituted in the free energy, 
the inverse capacitance can be written in terms of the G 
matrices 



— (-TOO — <-JOA(j^B(-rBO 



(27) 



The order of approximation in this method is the number 
of the highest multipoles 1,1' included. Matrix Gqo has 
dimension Ns x Ng with Ns the total number of conduc- 
tors. Matrices Gqa and Gbo are Ng x {NsNtotaimuiti) and 
{NsNtotahnuiti)xNs respectively, and matrix Gab has di- 
mension (NsNtotalmulti) X {N sNtotalmulti) ■ Ntotalmulti is 

the maximum number of multipolar terms considered. 
Formally it is 



^totalmulti — 



(28) 



with Imax the order of the maximum multipole included 
in the approximation. However the symmetries of the 
problem can help us to reduce it as the Gab elements 
corresponding to certain mi can be seen to vanish by 
symmetry. In particular in the case of azimuthal sym- 
metry, considered in the text, only m = gives non-zero 
values and the number of terms included can be reduced 
to Ntotalmulti = Imax- Depending on the geometry of the 
conductors it can be convenient to use different number 
of Imax for different conductors. In particular in the case 
of an array of small islands sandwiched by two large elec- 
trodes, it is better to use a larger number of multipoles 
at the electrodes. Most of the cases presented here are 
done with Imax ^ 8. 

The expression for G^^ follows from the decomposition 
of l/|a— b— R|, with a, b and R three points in space and 
depends on the geometry of the conductors. For a =^ (3 



G 



al3 

Iimil2m2 



{h + h - mi - m2)!(/i + I2 — mi + 1712)1 



{h + mi)!(/i - mi)!(/2 + ^2)1(^2 - 7712)! 

(-l)''+"'^il+i2+mi-m2(a^/3 " X 



mQl ■ 



(24) 



Substituting (^1)) into (^5)) and comparing it with ([T]), 
the inverse capacitance matrix can be expressed as 



■yv 



l^m,l' ,m' ,a,f3 



l,m,l,m' Lm I'n 



(25) 



with Ii rii the irregular solid spherical harmonics 
1 



/im(r) 



47r 



'+1 V 2/ + 1 



(30) 



The sign of G'^jln-^i^m^ depends not only on I2 and m2, 
but also on the order a(3 or j3a through the dependence 

Ili+l2+mi-m2{xi3 —Xa)- 
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For the case of an sphere a with radius Ra , " , ^ 

^ ti ,7711 ,t2?Tl2 



2/1- 



(31) 



The case of spheres on a row is especially simple. In 
this case we have azimuthal symmetry what means that 
all terms involving m ^ should vanish. Thus at order 
Imax we have just Ntotaimuiu = Imax- This simphfication 
allows us to go to reasonably high orders. We can elim- 
inate the indexes mi^m2 from the matrix G. Together 
with the diagonal terms G"" calculated above, and using 
that 



4tt 
21 + 1 



Pi [cost 



(32) 



and Pi{l) = 1 and Pi{-1) 
greatly simplified. Thus 



(—1)' the equations are 



^hM- IA1\ y-^) + l ' > 

for a ^ p. Here is the distance between the centers 
of the spheres a and (i. The diagonal of Go a and Gao 



are zero and G' 



and 



Q/3 

OA 



G%. Note that 



•-^00 - 



1 



(34) 



(35) 



The zero-order approximation recovers our expectation 
for the case of far-apart spheres. The correction to the 
inverse capacitance due to the higher order multipoles is 
given by — GqaG^^Gbo- As spheres come closer, higher 
order terms become more and more important. This is 
reasonable taking into account that the interaction be- 
tween two multipolar charges Qf^^^ and Qf „j decays as 

' af3 



IX. APPENDIX II. INTERACTION BETWEEN 
CHARGES AND CAPACITANCE MATRICES 

When two conductors become closer together the mo- 
bile charges in their surfaces screen the interaction be- 
tween the charges in them, compared to the 1 /r Coulomb 
law which describes the interaction between two isolated 
point-like charges. Other metallic systems in the sur- 
rounding environment contribute to this screening. It is 
expected that the charging energy of an sphere, or the 
energy to create an electron-hole pair between two is- 
lands will depend on the distance between the particles 
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(a) Island inverse capacitance C^^^, nearest-neighbor 
and excitonic energy, all in units of 
(C"°')~^, at the center of a 50 nanoparticle array (without 
electrodes) as a function of the interisland separation. Solid- 
lines give the value obtained with a 1/r interaction between 
charges. The effect of screening is evident for d/r"' < 1 — 2. 
All the plotted quantities saturate to a finite value as d/r"' 
vanishes, (b) Decay of the interaction potential from the center 
of a 50 nanoparticle array (without electrodes) as a function 
of island position for difi'erent values of d/r"'. Inverse capac- 
itances are given in units of (C"')~^ Solid lines correspond 
to an unscreened 1/r law. The bump is clearly observed only 
for very small d/r**'. It is much weaker than the one found 
inS and is restricted to the first nearest neighbor. For large 
distances the 1/r law is approached from above. The effect 
of screening is negligible for d/r ~ 3. In (a) and (b) dashed 
and dotted lines are included as a guide to the eye. 



in the array. In this appendix we describe the modifica- 
tion of the interaction due to screening and compare it 
with other models and calculations available in the lit- 
erature, is calculated as described in Appendix I. 
Both methods described give the same value of G~^ to 
several digits when the accuracy used in the computation 
is good enough. Here, inverse capacitances are measured 
in units of C^^l, the value of the inverse capacitance of 
a nanoparticle C~^^ when it is isolated. We find that the 
effect of screening is essentially restricted to the inter- 
actions between charges which are in the same islands 
or in the closest neighboring islands and only when the 
particles are very close d/r*^' < 1 — 2. 
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Three important quantities for the transport are the 
island inverse capacitance C^^ , the nearest-neighbor in- 
teraction C~Ij^i and the resulting excitonic energy E!l^'\ 
In Fig. 8(a) we show how these quantities depend on the 
distance between the particles, d/r"^"^ for the case of an 
array with N = 100 equal-sized islands, at the center of 
the array. The effects of screening start to be relevant 
for d/r'**' < 1 — 2. At this value C^^^ starts decreasing, 
and when d/r*'*' 0, approximately at d/r^^^ = 0.002 it 
saturates at about 0.68 of its bare value. C~Ij^i increases 
following a 1/r law as d/r is reduced, until d/r'^"^ 1 — 2. 
Then it decreases slightly due to the screening of the in- 
teraction, and finally it saturates at about d/r*''' =0.1. 
At large d/r*'*' the E^~^ starts decreasing when djr"^^^ 
is reduced, according to the 1/r increase of C~j^^. For 
fj^jj-isi ^ 2^ its value is affected both by the increase of 
C^i+i decrease of C,^^, at small distance is con- 

trolled by this last effect. Finally it saturates. 

The modification of the distance dependent interaction 
in one-dimensional arrays has been previously studied for 
some particular island geometries. For the case of thin 
circular disks (with disk axis perpendicular to the array 
axis). Whan et al.^* found that a 1/r law describes well 
the dependence of the non-diagonal inverse capacitance 
matrix elements on distance and that the diagonal ele- 
ments are reduced. Likharev and Matsuoka^'' analyzed 
the cases of an array of cubic islands (with inverse ca- 
pacitance matrix calculated using the finite-differences 
software FASTCAP) and a continuum model in which 
the discrete periodic structure is replaced by a continu- 
ous dielectric medium. They found that the dependence 
of the interaction on the distance between the charges 
r has a crossover between a linear dependence at short 
distances and a 1/r law at long distances. There is a 
bump on the interaction when compared to a 1/r law. 
The potential approaches the 1/r law from above. The 
interaction potential could by approximated by the ex- 
pression 



U{m) 



a 



a \mo 



-exp 



mo 



1 — exp 



mo 



(36) 

Here m is the distance in units of the array period, a, 
nT-o = To /a, with rp = Se/n a characteristic decay length 
of the interaction, e, the inter-islands dielectric constant, 
S the junction surface and n a fitting parameter with 
value very close to 1 and related to a by 



2 K 
K 2 



(37) 



The dependence of C~l^j on j that we have obtained 
for the experimentally relevant case of an array of spher- 
ical particles is plotted in Fig. 8(b). Screening is less 
important for the case of spheres than for cubic islands 
but larger than for the thin disks analyzed by Whan et 
al^. Compared with a 1/r law, at short distance the 
screened interaction potential decreases and at large dis- 
tances it increases. When the particles are very close. 



there is a bump in the renormalization of the interac- 
tion. Only when two charges are in the same particle 
or at the nearest neighbor one the interaction between 
them decreases. In other cases the interaction increase. 
We have not been able to fit our results with ((36l) . We re- 
cover the bump found by Likharev and Matsuoka— , but 
it is smaller than the one they observe, even when our 
particles are very close. We believe that the differences 
are due to the different geometry of the system under 
study. 

From the decomposition of the induced screening 
charge in high-order multipoles discussed in previous ap- 
pendix it is possible to get some insight on how does 
the bump appear. The correction to the inverse ca- 
pacitance SC~^ due to screening is —GoAG^gGg^, see 
((27| . which has the same sign as the monopole-induced 
multipole interaction. Thus, SC~} has the same sign 
as the interaction of the monopolar charge in a with 
the multipolar charge induced in all conductors by the 
monopolar charge in (3. Let us restrict to dipolar or- 
der, which will be the largest contribution. Qfj generates 
dipolar charges in all the other conductors 7, equal to 
QiP _ ^^{GJi )~^GioQi3 and the monopolar charge in 
a interact with these induced charges via Gqi- Both 
the monopolar-dipolar interaction Goi and the dipolar 
induced charges are odd quantities with respect to po- 
sition. To dipolar order, all the elements of G^^ are 
positive. The sign of the two odd quantities will control 
the sign of SG~p. 

Consider SG^^, with < a < (3 < N, the charges in- 
duced by a monopolar charge placed in /3 in conductors 
from 1 to a — 1 and from a + 1 to /? — 1 have the same 
sign, but the interaction of a monopolar charge in a with 
them, have opposite. Thus those terms coming from the 
charge induced in conductors from 1 to a — 1 and those 
from conductors a + 1 to f3 — 1 contribute to SC^^ with 
different sign. The sign of the contribution of the con- 
ductors at the right of (3 will be the same as those at the 
left of a and both induced charge and interaction have 
opposite sign. As further are these conductors to a and 
f3 the correction will be smaller. Individual contributions 
from each conductor 7 will be larger when a < 7 < /3. 
The contribution of the interaction term which comes 
from the dipoles generated in conductors between a and 
/3 increase G~^ . Only when the contribution of the terms 
which decrease the inverse capacitance matrix element is 
able to compensate the contribution of those ones which 
increase it, the total change will be negative. As closer 
are a and (3 the number of terms increasing the interac- 
tion will decrease. G~^ is expected to be smaller than the 
bare value, only if a and P are very close, what results in 
the appearance of the bump and the anti-screening effect 
at intermediate and large distances. 

The interaction obtained differs considerably from the 
one resulting from capacitive coupling only to nearest- 
neighbors. This form of the interaction has been used 
frequently in the literature, and applies when the system 
is coupled to a gate electrode which screens the long- 
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range part. The capacitance matrix has a triband form. 
Only Cii — Co + 2C and Ci^i±i = —C are finite. The 
inverse of Cy for the case of an infinite array becomes 

C^' = (38) 
with i and C°° defined from C°° = 2Cs/i(C"i) = 

1 /2 

(Cq + 4CCo) . Due to the exponential dependence ^ 
can be viewed as the decay length of the interaction, 
which increases with C/Cq. Interactions on the same 
island are given by C~j^ = 1/C°°. In the case of a finite 
array of length N this value is approached, from below, 
as N increases. The onsite case is recovered when (7 = 
and long-range interactions appear in the opposite limit 
Co/C ^ 0. In the later limit the interaction potential 
goes like 

C- = ^(f-|-.l) (39) 

It decays linearly with distance. The energy to create an 
electron-hole pair i^r" = l/2Cri + l/2C-±\,,±, - C-^i 
remains bound and equal to 1/ (2C) . On the contrary, 
the diagonal element C~j^ — diverges with the array 
size. There is not such unphysical divergence in C^^ with 
the array size in our model. 

To analyze the transport properties the array is sand- 
wiched between two electrodes, much larger than each 
of the nanoparticles. To this end we consider a one- 
dimensional array of N nanoparticles placed in between 
two large spheres, with radius R, which play the role 
of the leads. The large size ensures large screening and 
that Cqq and C^]^-^ are much smaller than the is- 
lands C~^^. The spherical shape greatly simplifies the 



calculations of the inverse capacitance matrix. The in- 
teraction between the charges at the islands and those at 
the electrodes and the inverse capacitance elements of the 
electrodes determine Af and Af, which control the po- 
larization voltage drop through the array and to a large 
extent the current flow at small voltages, see Figs. 6(a) 
and (b). 

For the size of the electrodes used in the text, R ~ 
50 — lOOr*'*', the inverse capacitance of the islands close 
to the electrodes is slightly reduced compared to those at 
the center, except for very small d/r*^'. For small d/r^"^ 
the inverse capacitance of islands at the center of the ar- 
ray is almost insensitive to the presence of the electrodes. 
If the electrodes are much larger the dependence of the 
inverse capacitance matrix elements with the size of the 
electrodes can become non-monotonous. This behavior, 
like the one found for small d/r"' is most probably as- 
sociated to the spherical shape chosen to model the elec- 
trodes. 
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